function LTP_Auto()
% Analyze the time course of EPSP changes
% Analysis includes: amplitude, slope, apparent IPSP, Ihold, and Vrmp.
% 11/1/99. Updates to clean up.
% LTP_Auto is taken from EPSP_TC, 2/28/03
% Does automated analysis, based on the filename of the selected data set.
%
global CONTROL

%try
    sf = getmainselection;

    if(sf > 0)
        pflag = getplotflag;
        QueMessage('LTP-Auto analysis', 1); % clear the que
        for i = 1:length(sf)
            sel = find_protocol(sf(i), 'timing_base'); % find protocol in this file
            if(~isempty(sel))
                fprintf(1, 'File: %s   index: %d\n', CONTROL(sel(1)).filename, sel(1));
                LTP_Auto2(sel(1), pflag);
            else
                fprintf(1, 'File: %s has no timing_base protocol\n', CONTROL(sf(i)).filename);
            end;

        end;
    end;
%catch
%    watchoff;
%    QueMessage('Error in LTP_Auto Analysis routine', 1);
%end

    function LTP_Auto2(sf, plot_flag)

        global VOLTAGE CURRENT
        global CONTROL DFILE
        % first generate empty output arrays

        err = 0;
        EPSPresult=[];

        QueMessage('LTP_Auto - Starting');
        dat = [];
        time = [];
        ivresult = []; % initialize it
        do_spike = 0;	% 0 turns off detailed spike analysis; 1 turns it on
        decay_start = 0.85; % position for decay fitting...
        fit_re_peak = 1; 	% do exp fits from peak of trace rather than "0" time

        [DFILE, err] = analysis_setup(DFILE, sf);

        if(err ~= 0)
            return;
        end;

        %-----------------------------Begin analysis--------------------------------

        % abstract general information
        spike_thresh=number_arg(CONTROL(sf).thresh); % to delete traces with spikes...
        QueMessage(sprintf('LTP_Auto: Spike Threshold: %7.2f mV', spike_thresh));
        protocol=deblank(lower(CONTROL(sf).protocol));
        rate = DFILE.rate.*DFILE.nr_channel/1000;
        [records,pts]=size(CURRENT);

        [RL, err] = record_parse(CONTROL(sf).reclist);
        frec=min(RL);
        lrec=max(RL);

        QueMessage('LTP_Auto - time base');
        % compute the time base for plotting (time is [rec, npts] array for EACH record)
        time=make_time(DFILE);
        tmax=max(max(time));
        k = find(DFILE.rate < 1);
        DFILE.rate(k) = 50;
        RATES = (DFILE.rate .* DFILE.nr_channel)/ 1000; % array of sampling rates, convert to msec


        % always compute the time of each stimulus, in seconds from the start of the data
        if(DFILE.mode >= 5)
            wz=DFILE.ztime;
            w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
            if(~isempty(w))
                wz(w+1:end)=wz(w+1:end)+24*60*60;
            end;
            zt = (wz-wz(1))/(60);
            cond_baseline = 5; % 5 min baseline
        else
            zt = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
            cond_baseline = 3;
        end;
        TM=zt;

        if(DFILE.mode >= 5)
            wz=DFILE.ztime;
            w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
            if(~isempty(w))
                wz(w+1:end)=wz(w+1:end)+24*60*60;
            end;
            ZT = (wz-wz(1))/(60);
            cond_baseline = 5; % 5 min baseline
        elseif length(DFILE.ztime) == length(RL)
            ZT = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
            cond_baseline = 3;
        else
            ZT = ones(length(RL), 1)*DFILE.cycle;
            cond_baseline = 5;
        end;

        % access ZT with the record number (1..n) to get the corresponding ztime
        % have to get the records first

        % Now,get the times when the valves switched (if any...)
        % and generate periods with valve 1, 2, 3 or 4.
        QueMessage('LTP-Auto - Valves...');
        p=datac('getnote'); % read the current notefile information.
        t_sw_valve=[];
        n_valve=[];
        TL=[];
        VL=[];
        if(~isempty(p) & length([p.proto]) > 0) % there should be some, but if not, don;'t do much
            % first set of arrays are immediate representations.
            sw_valve=[1 [p(find(diff([p.valve])~=0)+1).frec]];	% valve switch list (records)
            if(length(sw_valve) > 1)
                n_valve=[1 p(find(diff([p.valve])~=0)+1).valve]; % which valve...
                t_sw_valve=ZT(sw_valve); % don't forget offset from start of data..
                % now make long time arrays to match the other arrays.
                for i = frec:lrec
                    TL(i)=ZT(i);
                    for j=1:length(sw_valve)
                        if(i >= sw_valve(j))
                            VL(i)=n_valve(j);
                        end;
                    end;
                end;
            end;
        end;

        % Get analysis windows
        QueMessage('LTP-Auto - analysis windows');
        stim_list=CONTROL(sf).stim_time; % get the array
        psp_time=CONTROL(sf).psp_time; % get psp definition array
        if(ischar(psp_time))
            psp_time = str2num(psp_time);
        end;
        l = length(psp_time);
        if(l > 4)
            psp_time = reshape(psp_time, l/2, 2)';
        end;
        p2=size(psp_time);
        pair=(p2(1)>1);

        ts_rmp=number_arg(psp_time(1,1)); % get the baseline dur
        ts_b1=number_arg(psp_time(1,2)); % get start of psp list
        ts_e1=number_arg(psp_time(1,3)); % get end of psp list
        ts_i1=number_arg(psp_time(1,4)); % use s4 to find min voltage (IPSP max) after EPSP
        if(ts_e1 == ts_i1)
            ts_i1 = ts_e1+1;
        end;

        if(pair)
            ts_rmp2=number_arg(psp_time(2,1)); % get the baseline dur
            ts_b2=number_arg(psp_time(2,2)); % get start of psp list
            ts_e2=number_arg(psp_time(2,3)); % get end of psp list
            ts_i2=number_arg(psp_time(2,4)); % use s4 to find min
            if(ts_e2 == ts_i2)
                ts_i2 = ts_e2+1;
            end;
        end;

        % calculate measurement times
        QueMessage('EPSP_TC - measurement times');
        for i=1:records
            trmp(i) = floor(ts_rmp/RATES(i) - 0.5); % for RMP/ihold determination
            epsp_b1(i)=floor(ts_b1/RATES(i)-0.5); % window for begining of epsp
            epsp_e1(i)=floor(ts_e1/RATES(i)-0.5); % end of epsp
            ipsp_i1(i)=floor(ts_i1/RATES(i)-0.5); % end of ipsp that follows epsp.
            tmaxw(i)=ipsp_i1(i);	% max time in analysis

            % eventually we should make tsxx an array to generalize processing
            % but for first try, use fixed variables to get algorithm working.
            if(pair)
                epsp_b2(i)=floor(ts_b2/RATES(i)-0.5); % for peak iv
                epsp_e2(i)=floor(ts_e2/RATES(i)-0.5);
                ipsp_i2(i)=floor(ts_i2/RATES(i)-0.5);
                tmaxw(i)=ipsp_i2(i);
            end;
        end;

        % get step windows  (current pulse)
        ts1=number_arg(CONTROL(sf).durho);
        ts2=number_arg(CONTROL(sf).durs1);
        ts3=number_arg(CONTROL(sf).durs2);
        ts4=number_arg(CONTROL(sf).durs3);
        ts1_p = floor(ts1/RATES(1)-0.5); % pointer to end of durho
        ts2_p = floor((ts1+ts2)/RATES(1)-0.5); %pointer to end of s2
        ts3_p = floor((ts1+ts2+ts3)/RATES(1)-0.5); % pointer to end of s3
        ts4_p = floor((ts1+ts2+ts3+ts4)/RATES(1)-0.5); % pointer to end of s4
        if(ts1_p > pts)
            ts1_p = pts;
        end;
        if(ts2_p > pts)
            ts2_p = pts;
        end;
        if(ts3_p > pts)
            ts3_p = pts;
        end;
        if(ts4_p > pts)
            ts4_p = pts;
        end;



        %-----------------------------Begin analysis in earnest --------------------------------
        QueMessage('EPSP_TC - Measuring Ih and rmp');

        % measure holding current and "rmp"
        % also measure Rin with s1 pulse after ho (ts2)
        % and tau by fitting to s1 pulse.
        for i = 1:records
            hold_cur(i) = mean(CURRENT(i,1:trmp(i))); % save in array for later usage
            RMP(i) = mean(VOLTAGE(i,1:trmp(i)));
        end;
        CONTROL(sf).iHold = mean(hold_cur);
        CONTROL(sf).Rmp = mean(RMP);
        stim1 = CONTROL(sf).stim_time(1);

        % Now smooth the voltage out a bit
        QueMessage('EPSP_TC - Smoothing');
        for i = 1:records
            fsamp = 1000/RATES(i); % get sampling frequency
            fco = 1000;		% cutoff frequency in Hz
            wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
            if(wco < 1) % if wco is > 1 then this is not a filter!
                [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
                vsmo(i,:) = DigitalFilt(b, a, VOLTAGE(i,:)')'; % filter all the traces...
            else
                vsmo(i,:) = VOLTAGE(i,:);
            end
        end

        %vsmo = voltage; % in case you don't want filtering
        QueMessage('EPSP_TC - finding spikes');

        % find spikes in the records. We can't analyze an EPSP if spikes are present.
        [rmp_spike, first_isi, nr_spikes, spike_train] = find_spikes(DFILE, ts_b1, ts_e1, spike_thresh);
        srmp_list = find(isnan(rmp_spike) | rmp_spike >= ts_rmp); % NO spikes in rmp window
        srmp_nlist = find(~isnan(rmp_spike) | rmp_spike < ts_rmp); % spikes IN rmp window

        s1_list = []; s1_nlist = []; s2_list = []; s2_nlist = [];
        TM_SPIKE1=[]; TM_SPIKE2=[];

        for i = 1:records
            swin = floor([1:trmp(i) epsp_b1(i):ipsp_i1(i)]);
            if(any(VOLTAGE(i,swin)>spike_thresh))
                s1_list = [s1_list i];
            else
                s1_nlist = [s1_nlist i];
            end;
        end;
        TM_SPIKE1=TM(s1_nlist);

        if(pair)
            for i = 1:records
                swin = floor([1:trmp(i) epsp_b2(i):ipsp_i2(i)]);
                swin = floor([1:trmp(i) epsp_b2(i):ipsp_i2(i)]);
                if(any(VOLTAGE(i,swin)>spike_thresh))
                    s2_list = [s2_list i];
                else
                    s2_nlist = [s2_nlist i];
                end;
            end;
            TM_SPIKE2=TM(s2_nlist);
        else
            s2_list=1:records;
        end;

        QueMessage('EPSP_TC - Finding EPSP peak voltage and max rising slope');
        rtdur = floor((ts2_p-ts1_p)*0.8); % duration of window to use for rin measurement
        tr1 = ts1_p - rtdur; % baseline for voltage before step
        tr2 = ts2_p-rtdur;  % start of measurement in the pulse
        % get baseline voltage and ihold current
        for i = 1:records
            vbase(i)= mean(vsmo(i,1:trmp(i))); % mean(vsmo(i, (epsp_b1(i):epsp_b1(i)+3)))
            if(pair)
                vbase2(i)=mean(vsmo(i, (epsp_b2(i):epsp_b2(i)+3))); % mean(vsmo(i,1:trmp(i)));
            end;
            ibase(i)=mean(CURRENT(i,1:trmp(i)));

            if(ts1_p ~= ts2_p & rtdur > 0) % make sure valid values are available
                Rin(i) = - mean(VOLTAGE(i, tr2:ts2_p)) + mean(VOLTAGE(i, tr1:ts1_p)); % voltage difference
                rini =  - mean(CURRENT(i, tr2:ts2_p)) + mean(CURRENT(i, tr1:ts1_p));
                Rin(i) = Rin(i)/rini;

                [a1 a2 taum] = expfit(time(i, ts1_p+3:ts2_p)-time(i, ts1_p), VOLTAGE(i, ts1_p+3:ts2_p)); % time constant meausurement
                Tau(i) = taum;
            else
                Rin(i) = NaN;
                Tau(i) = NaN;
            end;
        end;
        Tau(find(Tau <= 1)) = NaN; % eliminate too-small tau's
        Tau(find(Tau > 5*ts2)) = NaN; % eliminate too-large taus

        % find EPSP amplitudes
        cmd = 2;
        rise = 4;
        decay = 10;
        iminus = rise*2;
        iplus = decay*5;
        SIGN=1;
        N = size(CURRENT);
        ihist = []; % peak current histogram data
        qhist = []; % total charge histogram data
        isihist = []; % inter-event ('spike') histogram data
        itau = [];
        itaum = [];
        iamp = [];
        icur = zeros(1,10); % current in acceptable traces
        m = 0; % count acceptable traces
        twin = [0 0];
        im = [0 0];
        stim_time = 5;

        k = 0;
        for i = s1_nlist
            k = k+ 1;
            [e_cand, e_amp, v_off, crit] = ClementsBekkers(VOLTAGE(k,:)-vbase(k), RATES(i), rise, decay, 0, 1, 0, 3000, 2);
            jl = 0; % index to last real event
            evcount = 0; % total count of events in this trace
            p=(DFILE.rate(k)*DFILE.nr_channel(k)/1000);
            npts=DFILE.nr_points;
            istd = std(VOLTAGE(k,[1:trmp-1])); % mean_std('current', k);
            istdn = -2*istd;
            istdp = 2*istd;
            lj = 0;
            QueMessage(sprintf('LTP_Auto analysis: %d candidate events', length(e_cand)));
            for i = 1:length(e_cand)-1 % go through all the candidate points
                j = e_cand(i);

                hold off;
                %      fprintf(1, '%d: j=%d jl=%d icand(i+1)=%d\n', i, j, jl, icand(i+1));
                if (((j - jl) > floor(iminus/p))  & (abs((e_cand(i+1)-j)) > floor(iplus/p) ))
                    prewin = j-[5:(floor(iminus/p))];
                    twin = ([-iminus:p:iplus]); % generates a time base for this plot
                    iwin = [1:length(twin)]+(j-floor(iminus/p));
                    % criteria for acceptance:
                    %    must fit in the data window completely,
                    %    there must be no large events greater than maxamp
                    if(min(prewin) > 0 & ...
                            min(iwin) > 0 & ...
                            max(iwin) < npts) % & ...
                        %               all(abs(CURRENT(k,iwin)-imean) < par.maxamp)); % & ...
                        twin = ([-iminus:p:iplus]); % generates a time base for this plot
                        qwin = j + [-1:floor(2/p)];
                        iwin = [1:length(twin)]+(j-floor(iminus/p));
                        bl=mean(VOLTAGE(k,prewin));
                        iplot = VOLTAGE(k,iwin)-bl; % remove the baseline once and for all
                        go = 1; % final inclusion criterion - no junk going the other way
                        if((SIGN == 1) & any(iplot <= -istdn))
                            %     go = 0;
                        end;
                        if((SIGN == 0) & any(iplot >= istdp))
                            %   go = 0;
                        end;
                        if(go == 1)
                            if(SIGN == 1) % looking at positive-going events
                                [imax, jmax] = max(iplot);
                                ihist = [ihist imax];
                                qhist = [qhist sum(VOLTAGE(k,qwin)-bl)];
                                iamp = ihist;
                                i95 = decay_start*imax;
                                jmx = floor((length(iplot) - jmax)/2) + jmax;
                                [i95f p95f] = findnear(iplot(jmax:jmx), i95); % find point 95% way down falling slope
                                p95f = p95f(2)+jmax - 1;
                                if(fit_re_peak)
                                    tpk = twin(jmax);
                                else
                                    tpk = 0;
                                end;

                                if(jmx-p95f < 2)
                                    tau = NaN;
                                else
                                    [a0 a1 tau] = expfit(twin(p95f:jmx)-tpk, iplot(p95f:jmx)); % chebyshev fit.....
                                    [fpar, chisq, niter, volt_fit1] = mrqfit('exponential', [a0 a1 tau], twin(p95f:jmx)-tpk, iplot(p95f:jmx), [], ...
                                        [], [], [], 50, []);
                                    if(isempty(volt_fit1))
                                        fiti = a0 + a1*exp(-time_fit/tau);
                                        fitmeth = 'c'; % chebyshev
                                    else
                                        a0 = fpar(1); a1=fpar(2); tau = fpar(3);
                                        fiti = a1*exp(-(twin(p95f:jmx)-tpk)/tau)+a0;
                                        fitmeth = 'm'; % marquardt
                                    end;

                                end;
                                if(tau < p | tau > 2*iplus) tau = NaN; end; % limit taus to reasonable/acceptable values
                                itau = [itau tau];
                                itaum = [itaum fitmeth];
                                % fiti = a1*exp(-tn(p95f:end)/tau)+a0;
                            else % looking at negative going events
                                [imin, jmin] = min(iplot);
                                ihist = [ihist imin];
                                iamp = ihist;
                                qhist = [qhist sum(VOLTAGE(k,qwin)-bl)];
                                i95 = decay_start*imin;
                                jmx = floor((length(iplot) - jmin)/2) + jmin;
                                [i95f p95f] = findnear(iplot(jmin:jmx), i95); % find point 95% way down falling slope
                                p95f = p95f(2)+jmin - 1;
                                if(fit_re_peak)
                                    tpk = twin(jmin);
                                else
                                    tpk = 0;
                                end;
                                if(jmx-p95f < 2)
                                    tau = NaN;
                                else
                                    [a0 a1 tau] = expfit(twin(p95f:jmx)-tpk, iplot(p95f:jmx)); % chebyshev fit.....
                                    [fpar, chisq, niter, volt_fit1] = mrqfit('exponential', [a0 a1 tau], twin(p95f:jmx)-tpk, iplot(p95f:jmx), [], ...
                                        [], [], [], 50, []);
                                    if(isempty(volt_fit1))
                                        fiti = a0 + a1*exp(-time_fit/tau);
                                        fitmeth = 'c';
                                    else
                                        a0 = fpar(1); a1=fpar(2); tau = fpar(3);
                                        fiti = a1*exp(-(twin(p95f:jmx)-tpk)/tau)+a0;
                                        fitmeth = 'm';
                                    end;
                                end;
                                if(tau < p| tau > 2*iplus) tau = NaN; end; % limit taus to reasonable/acceptable values
                                itau = [itau tau];
                                itaum = [itaum fitmeth];
                            end;

                            h = findobj('Tag', 'Trace_analysis'); % check for pre-existing window
                            if(isempty(h))
                                h=figure('Tag', 'Trace_analysis');
                            end;
                            figure(h);
                            line(twin, iplot);
                            m = m + 1;
                            if(m == 1)
                                icur = zeros(1,length(iplot));
                            end;
                            icur(m,:) = iplot;
                            if(cmd == 2)
                                isamp(m) = e_amp(i);
                                icoff(m) = v_off(i);
                                ibl(m) = bl;
                            end;
                            hold on
                        end;
                    end;
                end;
                evcount = evcount + 1; % keep track of events (even if they are not included in traces...)
                if(evcount > 1)
                    isihist = [isihist time(j)-time(lj)]; % interevent distribution
                end;
                lj = j; % keep running track of most recent event
                jl = j+floor(iplus/p);
            end;

        end;


        [V_epsp_1(k) I_epsp_1(k)] = max(vsmo(i,epsp_b1(i):ipsp_i1(i))'-vbase(i));
        V_epsp_1(k) = mean(vsmo(i,(I_epsp_1(k)-2+epsp_b1(i)):(I_epsp_1(k)+2+epsp_b1(i)))'-vbase(i)); % smooth it out w/ 5 points
        %	V_epsp_1(k) = V_epsp_1(k) - vbase(i);
        T_epsp_1(k) = TM(i); % make matching time base

        slx = 1;
        [S_epsp_1(k)] = mean(diff(vsmo(i,epsp_b1(i):epsp_e1(i))))/RATES(i); % slope
        if(~isempty(DFILE.shkl))
            O_epsp_1(k) = (slx+epsp_b1(i))*RATES(i)-DFILE.shkl; % save onset time.
        else
            O_epsp_1(k) = (slx+epsp_b1(i))*RATES(i)-stim1; % save onset time.
        end;
        %S_epsp_1(k) = mean(diff(vsmo(i,slx-2+epsp_b1(i):slx+2+epsp_b1(i))'))/RATES(i); % slope
        T_s_epsp_1(k) = TM(i); % timing for slopes...

        [V_ipsp_1(k) I_ipsp_1(k)] = min(vsmo(i, epsp_e1(i):ipsp_i1(i))'); % ipsps.
        %	V_ipsp_1(k) = V_ipsp_1(k) - vbase(i);
        V_ipsp_1(k) = -mean((vsmo(i,(I_ipsp_1(k)-2+epsp_e1(i)):(I_ipsp_1(k)+2+epsp_e1(i)))'-vbase(i))); % smooth it out w/ 5 points
        T_ipsp_1(k) = TM(i);

        [S_ipsp_1(k) slx] = min(diff(vsmo(i, epsp_e1(i):(I_ipsp_1(k)+epsp_e1(i)))')); % slope - find point of max slope on RISING phase of IPSP
        S_ipsp_1(k) = -mean(diff(vsmo(i,slx-2+epsp_e1(i):slx+2+epsp_e1(i))'))/RATES(i); % slope - average around the pont
        T_s_ipsp_1(k) = TM(i); % timing for slopes...

  %  end;
    V_epsp_1(find(V_epsp_1<0)) = NaN; % clear out useless information
    S_epsp_1(isnan(V_epsp_1)) = NaN;
    S_epsp_1(find(S_epsp_1<0)) = NaN;
    O_epsp_1(isnan(S_epsp_1)) = NaN;

    V_ipsp_1(find(V_ipsp_1<0)) = NaN;
    S_ipsp_1(isnan(V_ipsp_1)) = NaN;
    S_ipsp_1(find(S_ipsp_1<0)) = NaN;

    if(pair) % second psp etc...
        k = 0;
        for i = s2_nlist
            k = k + 1;
            [V_epsp_2(k) I_epsp_2(k)] = max(vsmo(i,epsp_b2(i):epsp_e2(i))');
            %		V_epsp_2(k) = V_epsp_2(k) - vbase2(i);
            V_epsp_2(k) = mean(abs(vsmo(i,I_epsp_2(k)-2+epsp_b2(i):I_epsp_2(k)+2+epsp_b2(i))'-vbase(i))); % smooth it out w/ 5 points
            T_epsp_2(k) = TM(i); % make matching time base
            S_epsp_2(k) = max(diff(vsmo(i,epsp_b2(i):epsp_e2(i))')); % slope
            T_s_epsp_2(k) = TM(i); % timing for slopes...

        end;
    end;


    TMBL=TM(srmp_list);
    VBL=vbase(srmp_list); % our baselines are also abstracted this way...
    IBL=ibase(srmp_list);
    RINL = Rin(srmp_list);
    TAUL = Tau(srmp_list);

    % determine records for EPSP averaging and measurement, etc.
    % basic window is 3 - 0 min prior to first valve change, and last 3 min's
    % of the associated shock protocol for that valve change.
    % we use the exclusion lists to actually generate a short list to use.

    sw1 = []; sw2 = [];
    rv1 = []; tv1 = []; iv0 = []; rv0 = []; ep1l = [];
    if(~isempty(n_valve))
        sw1=sw_valve(find(n_valve==1)); % records when valve switched to 1.
        sw2=sw_valve(find(n_valve==2)); % records when valve switched to 2.
    end;
    if(~isempty(sw1) & ~isempty(sw2))
        ep1 = [sw1(1):sw2(1)-1];
        if(length(sw1) > 1)
            ep2 = [sw2(1):sw1(2)-1];
        else
            ep2 = sw2(1);
        end;

        rv1=sw_valve(find(diff(n_valve)==1)+1); % record and time at which first occured
    else
        sw1=0; sw2 = 0; ep1 = 0; ep2=0; rv1=1;
    end;
    tv1=ZT(rv1(1)); % use only the first one...
    iv0=find(ZT>(tv1-3*60)); % convert from mins to secs, get the time
    tv0=ZT(iv0(1));
    rv0=iv0(1);
    if(sw1 == 0 & sw2 == 0)
        sw1 = 1;
        sw2 = max(s1_nlist);
    end;
    ep1l=find(s1_nlist >= sw1(1) & s1_nlist < sw2(1));



    k = 0;
    for i = s1_nlist(ep1l)
        k = k + 1;
        %	fprintf('[1] i: %d  maxv: %6.1f\n', i, max(vsmo(i,epsp_b1(i):epsp_e1(i))));
        EPSP1(k,:)=vsmo(i, 1:ipsp_i1(i)) - vbase(i);
        T_EPSP1(k,:)=time(i, 1:ipsp_i1(i));
    end;
    % rotate before plotting;
    EPSP1 = EPSP1';
    T_EPSP1 = T_EPSP1';

    if(pair)
        ep2l=find(s2_nlist >= sw1(1) & s2_nlist < sw2(1));
        k = 0;
        for i = s2_nlist(ep2l)
            k = k + 1;
            fprintf('[2] i: %d  maxv: %6.1f\n', i, max(vsmo(i,epsp_b2(i):epsp_e2(i))));
            EPSP2(k,:)=vsmo(i, epsp_b2(i):epsp_e2(i)) - vbase2(i);
            T_EPSP2(k,:)=time(i, epsp_b2(i):epsp_e2(i))-time(i,epsp_b2(i));
        end;
        % rotate before plotting;
        EPSP2 = EPSP2';
        T_EPSP2 = T_EPSP2';
    end;

    QueMessage('EPSP_TC - analysis complete');



%-----------------------------Prepare for plotting--------------------------------

% for plotting, do baselines/stddev.
bl = find(T_epsp_1 < 5);
mv_e1 = mean(V_epsp_1(bl));
std_e1 = std(V_epsp_1(bl));
pstd_e1 = mv_e1 + std_e1;
nstd_e1 = mv_e1 - std_e1;

ms_e1 = mean(S_epsp_1(bl));
stds_e1 = std(S_epsp_1(bl));
pstds_e1 = ms_e1 + stds_e1;
nstds_e1 = ms_e1 - stds_e1;
tplot = [T_epsp_1(1)+0.1,T_epsp_1(end)-0.1];

mv_i1 = mean(V_ipsp_1(bl));
std_i1 = std(V_ipsp_1(bl));
pstd_i1 = mv_i1 + std_i1;
nstd_i1 = mv_i1 - std_i1;

ms_i1 = mean(S_ipsp_1(bl));
stds_i1 = std(S_ipsp_1(bl));
pstds_i1 = ms_i1 + stds_i1;
nstds_i1 = ms_i1 - stds_i1;

% save the results
if(DFILE.mode >= 5)
    PSP_TC.zt = DFILE.ztime; % we need ztime to align multiple traces...
else
    PSP_TC.zt = DFILE.ztime / (1000);
end;

% store the results.
PSP_TC.T_epsp_1 = T_epsp_1;
PSP_TC.V_epsp_1 = V_epsp_1;
PSP_TC.T_s_epsp_1 = T_s_epsp_1;
PSP_TC.S_epsp_1 = S_epsp_1;
PSP_TC.T_ipsp_1 = T_ipsp_1;
PSP_TC.V_ipsp_1 = V_ipsp_1;
PSP_TC.O_epsp_1 = mean_var(O_epsp_1); % just return the mean max slope time.

% baseline measurements for comparison...
PSP_TC.T = TMBL;
PSP_TC.Vm = VBL;
PSP_TC.Rin = RINL;
PSP_TC.Ih = IBL;
PSP_TC.Tau = TAUL;

% and average EPSC during the last 2 minutes of the sample period.

CONTROL(sf).PSP_TC = PSP_TC;

%----------------------------- plot if figure is set--------------------------------


